import psycopg2
import psycopg2.extras
import sys
import gdal
from gdalconst import *
import os
import math
from shapely.wkt import dumps, loads
import datetime
from datetime import timedelta
import numpy as np

def interpolate(x1,y1,x0,y0,x):
    y=math.exp((x-x0)*((math.log(y1)-math.log(y0))/(x1-x0))+math.log(y0))
    return y

def corel_mod04_aeronet(time_range):
    in_file = "correl/mod04_aronet_2.csv"
    station_query= 'SELECT id, ST_AsText(location) as location FROM public.grdstation WHERE id between 68 and 74'
    mod_query="SELECT aqstime, st_neighborhood(rasref, (st_GeomFromText('POINT({0} {1})', 4326)), 2, 2, true) as aod FROM apom.satresampmod04 as mod04 order by mod04.aqstime asc"
    aeronet_query="SELECT avg(aod_675) as aod675,avg(aod_500) as aod500 from apom.grdaeraot_data where processtype=2 and (aod_500>=0 and aod_675>=0) and aqstime > '{0}'::timestamp - interval '{1} minutes' and aqstime < '{2}'::timestamp + interval '{3} minutes' and stationid = {4} GROUP BY aqstime"
    f = open(in_file, "w+")
    con = None
    try:
        # doc danh sach tram
        con = psycopg2.connect(host='192.168.0.4', database='fimo', user='rasdaman', password='rasdaman') 
        cur = con.cursor(cursor_factory=psycopg2.extras.DictCursor)
        #cur.execute('SELECT id, ST_AsText(location) as location FROM public.grdstation WHERE type = 2')
        cur.execute(station_query)
       
        stations = cur.fetchall()
        for station in stations:
            stationid = station["id"]                
            location = loads(station["location"])
            # lay cac pixel xung quanh tram tren anh mod04,07
       
            cur.execute(mod_query.format(location.x, location.y))
           
          
            rows = cur.fetchall()
            for row in rows:
                aqstime = datetime.datetime.strptime(str(row["aqstime"]), "%Y-%m-%d %H:%M:%S")
                #aqstime = datetime.datetime.strptime(str(row["aqstime"]), "%Y-%m-%d %H:%M:%S") + timedelta(hours=7)
    
                # du lieu gia ve la 1 mang 2 chieu
                aod = row["aod"]
                aod_result = []
                for i in range(len(aod)):
                    # filter nhung gia tri != None thi se them vao mang ket qua
                    tmp = filter(lambda a: a != None, aod[i])
                    aod_result += tmp
                # calculate real aod value
                mod04_aod = (np.median(aod_result) * 0.00100000004749745) if len(aod_result) > 0 else 0
                    
                if mod04_aod > 0:
                    
                    cur.execute(aeronet_query.format(aqstime, time_range, aqstime, time_range, stationid))                  
                    pm_row = cur.fetchone()
                    if not(pm_row is None):
                        aronet_aod=interpolate(675, pm_row["aod675"], 500, pm_row["aod500"], 550)
                        f.write("{0},{1},{2},{3},{4},{5}\n".format(stationid, aqstime.year, aqstime.month, aqstime, aronet_aod, mod04_aod))
                #print aqstime
    except psycopg2.DatabaseError, e:
        print 'Error %s' % e    
        sys.exit(1)
    finally:
        if con:
            con.close()
        if f:
            f.close()
    print "done"
    data = np.loadtxt("correl/mod04_aronet_2.csv", delimiter = ",", usecols=(4,5))
    a= np.corrcoef(data[:,0],data[:,1])
    print "number of sample",len(data)
    print "mod04-aeronet R2:",a[1,0]*a[1,0]
def corel_myd04_aeronet(time_range):
    in_file = "correl/myd04_aronet_2.csv"
    station_query= 'SELECT id, ST_AsText(location) as location FROM public.grdstation WHERE id between 68 and 74'
    mod_query="SELECT aqstime, st_neighborhood(rasref10, (st_GeomFromText('POINT({0} {1})', 4326)), 2, 2, true) as aod FROM apom.satresampmyd04 as mod04 order by mod04.aqstime asc"
    aeronet_query="SELECT avg(aod_675) as aod675,avg(aod_500) as aod500 from apom.grdaeraot_data where processtype=2 and (aod_500>=0 and aod_675>=0) and aqstime > '{0}'::timestamp - interval '{1} minutes' and aqstime < '{2}'::timestamp + interval '{3} minutes' and stationid = {4} GROUP BY aqstime"
    f = open(in_file, "w+")
    con = None
    try:
        # doc danh sach tram
        con = psycopg2.connect(host='192.168.0.4', database='fimo', user='rasdaman', password='rasdaman') 
        cur = con.cursor(cursor_factory=psycopg2.extras.DictCursor)
        #cur.execute('SELECT id, ST_AsText(location) as location FROM public.grdstation WHERE type = 2')
        cur.execute(station_query)
       
        stations = cur.fetchall()
        for station in stations:
            stationid = station["id"]                
            location = loads(station["location"])
            # lay cac pixel xung quanh tram tren anh mod04,07
       
            cur.execute(mod_query.format(location.x, location.y))
           
          
            rows = cur.fetchall()
            for row in rows:
                aqstime = datetime.datetime.strptime(str(row["aqstime"]), "%Y-%m-%d %H:%M:%S")
                #aqstime = datetime.datetime.strptime(str(row["aqstime"]), "%Y-%m-%d %H:%M:%S") + timedelta(hours=7)
    
                # du lieu gia ve la 1 mang 2 chieu
                aod = row["aod"]
                aod_result = []
                for i in range(len(aod)):
                    # filter nhung gia tri != None thi se them vao mang ket qua
                    tmp = filter(lambda a: a != None, aod[i])
                    aod_result += tmp
                # calculate real aod value
                mod04_aod = (np.median(aod_result) * 0.00100000004749745) if len(aod_result) > 0 else 0
                    
                if mod04_aod > 0:
                    
                    cur.execute(aeronet_query.format(aqstime, time_range, aqstime, time_range, stationid))                  
                    pm_row = cur.fetchone()
                    if not(pm_row is None):
                        aronet_aod=interpolate(675, pm_row["aod675"], 500, pm_row["aod500"], 550)
                        f.write("{0},{1},{2},{3},{4},{5}\n".format(stationid, aqstime.year, aqstime.month, aqstime, aronet_aod, mod04_aod))
                #print aqstime
    except psycopg2.DatabaseError, e:
        print 'Error %s' % e    
        sys.exit(1)
    finally:
        if con:
            con.close()
        if f:
            f.close()
    print "done"
    data = np.loadtxt("correl/myd04_aronet_2.csv", delimiter = ",", usecols=(4,5))
    a= np.corrcoef(data[:,0],data[:,1])
    print "number of sample",len(data)
    print "mod04-aeronet R2:",a[1,0]*a[1,0]
def corel_myd04_aeronet_daily(time_range):
    in_file = "correl/myd04_aronet.csv"
    station_query= 'SELECT id, ST_AsText(location) as location FROM public.grdstation WHERE id between 68 and 74'
    mod_query="SELECT aqstime, st_neighborhood(rasref10, (st_GeomFromText('POINT({0} {1})', 4326)), 2, 2, true) as aod FROM apom.satresampmyd04 as mod04 order by mod04.aqstime asc"
    aeronet_query="SELECT avg(aod_675) as aod675,avg(aod_500) as aod500 from apom.grdaeraot_data where processtype=2 and (aod_500>=0 and aod_675>=0) and date(aqstime) = '{0}' and stationid = {1} GROUP BY date(aqstime)"
    f = open(in_file, "w+")
    con = None
    try:
        # doc danh sach tram
        con = psycopg2.connect(host='192.168.3.4', database='fimo', user='rasdaman', password='rasdaman') 
        cur = con.cursor(cursor_factory=psycopg2.extras.DictCursor)
        #cur.execute('SELECT id, ST_AsText(location) as location FROM public.grdstation WHERE type = 2')
        cur.execute(station_query)
       
        stations = cur.fetchall()
        for station in stations:
            stationid = station["id"]                
            location = loads(station["location"])
            # lay cac pixel xung quanh tram tren anh mod04,07
       
            cur.execute(mod_query.format(location.x, location.y))
           
          
            rows = cur.fetchall()
            for row in rows:
                aqstime = datetime.datetime.strptime(str(row["aqstime"]), "%Y-%m-%d %H:%M:%S") + timedelta(hours=7)
    
                # du lieu gia ve la 1 mang 2 chieu
                aod = row["aod"]
                aod_result = []
                for i in range(len(aod)):
                    # filter nhung gia tri != None thi se them vao mang ket qua
                    tmp = filter(lambda a: a != None, aod[i])
                    aod_result += tmp
                # calculate real aod value
                mod04_aod = (np.median(aod_result) * 0.00100000004749745) if len(aod_result) > 0 else 0
                    
                if mod04_aod > 0:
                    
                    cur.execute(aeronet_query.format(aqstime, stationid))                  
                    pm_row = cur.fetchone()
                    if not(pm_row is None):
                        aronet_aod=interpolate(675, pm_row["aod675"], 500, pm_row["aod500"], 550)
                        f.write("{0},{1},{2},{3},{4},{5}\n".format(stationid, aqstime.year, aqstime.month, aqstime, aronet_aod, mod04_aod))
                #print aqstime
    except psycopg2.DatabaseError, e:
        print 'Error %s' % e    
        sys.exit(1)
    finally:
        if con:
            con.close()
        if f:
            f.close()
    print "done"
    data = np.loadtxt("correl/myd04_aronet_2.csv", delimiter = ",", usecols=(4,5))
    a= np.corrcoef(data[:,0],data[:,1])
    print "number of sample",len(data)
    print "mod04-aeronet R2:",a[1,0]*a[1,0]


def corel_mod04_pm(time_range):
    f = open("correl/mod04_pm.csv", "w+")
    station_query='SELECT id, ST_AsText(location) as location FROM public.grdstation WHERE id >=1 and id <=5'
    mod_query="SELECT mod04.aqstime, mod04.aod from (SELECT aqstime, st_neighborhood(rasref10, (st_GeomFromText('POINT({0} {1})', 4326)), 2, 2, true) as aod FROM apom.satresampmod04 where aqstime < '2014-09-10 00:00:00'::timestamp) as mod04 order by mod04.aqstime asc"
    pm_query="SELECT avg(pm1), avg(pm25), avg(pm10) FROM apom.grdpmhiscem_data WHERE (pm1 > 0 OR pm25 > 0 OR pm10 > 0) and aqstime > '{0}'::timestamp - interval '{1} minutes' and aqstime < '{2}'::timestamp + interval '{3} minutes' and stationid = {4} GROUP BY aqstime"
    con = None
    try:
        # doc danh sach tram
        con = psycopg2.connect(host='192.168.3.4', database='fimo', user='rasdaman', password='rasdaman') 
        cur = con.cursor(cursor_factory=psycopg2.extras.DictCursor)
        cur.execute(station_query)
        stations = cur.fetchall()
        for station in stations:
            stationid = station["id"]
            # bo qua tram Khanh Hoa
            # if stationid  == 3:
            #     continue
                
            location = loads(station["location"])
    
            cur.execute(mod_query.format(location.x, location.y))
         
            rows = cur.fetchall()
            for row in rows:
                aqstime = datetime.datetime.strptime(str(row["aqstime"]), "%Y-%m-%d %H:%M:%S") + timedelta(hours=7)
            
                # du lieu gia ve la 1 mang 2 chieu
                aod = row["aod"]
                aod_result = []
                for i in range(len(aod)):
                    # filter nhung gia tri != None thi se them vao mang ket qua
                    tmp = filter(lambda a: a != None, aod[i])
                    aod_result += tmp
                # calculate real aod value
                aod_value = (np.median(aod_result) * 0.00100000004749745) if len(aod_result) > 0 else 0
                if (aod_value > 0):
                    cur.execute(pm_query.format(aqstime, time_range, aqstime, time_range, stationid))
                    pm_row = cur.fetchone()
                    if not(pm_row is None):
                        pm1 = pm_row[0]
                        pm25 = pm_row[1]
                        pm10 = pm_row[2]                 
                        f.write("{0},{1},{2},{3},{4},{5},{6},{7}\n".format(stationid, aqstime.year, aqstime.month, aqstime, pm1, pm25, pm10, aod_value))
            
    except psycopg2.DatabaseError, e:
        print 'Error %s' % e    
        sys.exit(1)
    finally:
        if con:
            con.close()
        if f:
            f.close()
    data = np.loadtxt("correl/mod04_pm.csv", delimiter = ",", usecols=(4,5,6,7))
    print "Number of sample:",len(data)
    
    a= np.corrcoef([data[:,0],data[:,3]])
    print "MOD04 - PM1 R2:", a[1,0]*a[1,0]
     
    b= np.corrcoef([data[:,1],data[:,3]])
    print "MOD04 - PM25 R2:", b[1,0]*b[1,0]
     
    b= np.corrcoef([data[:,2],data[:,3]])
    print "MOD04 - PM10 R2:", b[1,0]*b[1,0]
def corel_myd04_pm(time_range):
    f = open("correl/myd04_pm.csv", "w+")
    station_query='SELECT id, ST_AsText(location) as location FROM public.grdstation WHERE id >=1 and id <=5'
    mod_query="SELECT mod04.aqstime, mod04.aod from (SELECT aqstime, st_neighborhood(rasref10, (st_GeomFromText('POINT({0} {1})', 4326)), 2, 2, true) as aod FROM apom.satresampmyd04 where aqstime < '2014-09-10 00:00:00'::timestamp) as mod04 order by mod04.aqstime asc"
    pm_query="SELECT avg(pm1), avg(pm25), avg(pm10) FROM apom.grdpmhiscem_data WHERE (pm1 > 0 OR pm25 > 0 OR pm10 > 0) and aqstime > '{0}'::timestamp - interval '{1} minutes' and aqstime < '{2}'::timestamp + interval '{3} minutes' and stationid = {4} GROUP BY aqstime"
    con = None
    try:
        # doc danh sach tram
        con = psycopg2.connect(host='192.168.3.4', database='fimo', user='rasdaman', password='rasdaman') 
        cur = con.cursor(cursor_factory=psycopg2.extras.DictCursor)
        cur.execute(station_query)
        stations = cur.fetchall()
        for station in stations:
            stationid = station["id"]
            # bo qua tram Khanh Hoa
            # if stationid  == 3:
            #     continue
                
            location = loads(station["location"])
    
            cur.execute(mod_query.format(location.x, location.y))
         
            rows = cur.fetchall()
            for row in rows:
                aqstime = datetime.datetime.strptime(str(row["aqstime"]), "%Y-%m-%d %H:%M:%S") + timedelta(hours=7)
            
                # du lieu gia ve la 1 mang 2 chieu
                aod = row["aod"]
                aod_result = []
                for i in range(len(aod)):
                    # filter nhung gia tri != None thi se them vao mang ket qua
                    tmp = filter(lambda a: a != None, aod[i])
                    aod_result += tmp
                # calculate real aod value
                aod_value = (np.median(aod_result) * 0.00100000004749745) if len(aod_result) > 0 else 0
                if (aod_value > 0):
                    cur.execute(pm_query.format(aqstime, time_range, aqstime, time_range, stationid))
                    pm_row = cur.fetchone()
                    if not(pm_row is None):
                        pm1 = pm_row[0]
                        pm25 = pm_row[1]
                        pm10 = pm_row[2]                 
                        f.write("{0},{1},{2},{3},{4},{5},{6},{7}\n".format(stationid, aqstime.year, aqstime.month, aqstime, pm1, pm25, pm10, aod_value))
            
    except psycopg2.DatabaseError, e:
        print 'Error %s' % e    
        sys.exit(1)
    finally:
        if con:
            con.close()
        if f:
            f.close()
    data = np.loadtxt("correl/myd04_pm.csv", delimiter = ",", usecols=(4,5,6,7))
    print "Number of sample:",len(data)
    
    a= np.corrcoef([data[:,0],data[:,3]])
    print "MYD04 - PM1 R2:", a[1,0]*a[1,0]
     
    b= np.corrcoef([data[:,1],data[:,3]])
    print "MYD04 - PM25 R2:", b[1,0]*b[1,0]
     
    b= np.corrcoef([data[:,2],data[:,3]])
    print "MYD04 - PM10 R2:", b[1,0]*b[1,0]
def corel_mod07press_pm(time_range):
    f = open("correl/mod07press_pm_2.csv", "w+")
    station_query = 'SELECT id, ST_AsText(location) as location FROM public.grdstation WHERE id >=1 and id <=5'
    mod_query="SELECT aqstime, ST_Value(rasref10, (st_GeomFromText('POINT({0} {1})', 4326)),true) as pressure FROM apom.satresampmod07pressure where aqstime < '2014-09-10 00:00:00'::timestamp"
    pm_query="SELECT avg(pm1), avg(pm25), avg(pm10),avg(barometer) FROM apom.grdpmhiscem_data WHERE (pm1 > 0 OR pm25 > 0 OR pm10 > 0) and aqstime > '{0}'::timestamp - interval '{1} minutes' and aqstime < '{2}'::timestamp + interval '{3} minutes' and stationid = {4} GROUP BY aqstime"
    con = None
    try:
        # doc danh sach tram
        con = psycopg2.connect(host='192.168.3.4', database='fimo', user='rasdaman', password='rasdaman') 
        cur = con.cursor(cursor_factory=psycopg2.extras.DictCursor)
        cur.execute(station_query)
        stations = cur.fetchall()
        for station in stations:
            stationid = station["id"]
            # bo qua tram Khanh Hoa
            # if stationid  == 3:
            #     continue
                
            location = loads(station["location"])
    
            cur.execute(mod_query.format(location.x, location.y))                                                                                                                                                                                                                  
            rows = cur.fetchall()
            
            for row in rows:
                press_aqstime = datetime.datetime.strptime(str(row["aqstime"]), "%Y-%m-%d %H:%M:%S") + timedelta(hours=7)
                pressure = row["pressure"]
                if(pressure!=None):
                    pressure_value = pressure*0.100000001490116
                else:
                    pressure_value=0;

       
                if pressure_value > 0:
                 
                    cur.execute(pm_query.format(press_aqstime, time_range, press_aqstime, time_range, stationid))
                    
                    pm_row = cur.fetchone()
                    if not(pm_row is None):
                        pm1 = pm_row[0]
                        pm25 = pm_row[1]
                        pm10 = pm_row[2]
                        barometer=pm_row[3]
                      
                        f.write("{0},{1},{2},{3},{4},{5},{6},{7},{8}\n".format(stationid, press_aqstime.year, press_aqstime.month, press_aqstime, pm1, pm25, pm10,barometer,pressure_value))
            
    except psycopg2.DatabaseError, e:
        print 'Error %s' % e    
        sys.exit(1)
    finally:
        if con:
            con.close()
        if f:
            f.close()
    data = np.loadtxt("correl/mod07press_pm.csv", delimiter = ",", usecols=(4,5,6,7,8))
    print "Number of sample:",len(data)
    
    a= np.corrcoef([data[:,0],data[:,4]])
    print "MOD07press - PM1 R2:", a[1,0]*a[1,0]
     
    b= np.corrcoef([data[:,1],data[:,4]])
    print "MOD07press - PM25 R2:", b[1,0]*b[1,0]
     
    b= np.corrcoef([data[:,2],data[:,4]])
    print "MOD07press - PM10 R2:", b[1,0]*b[1,0]
    
    b= np.corrcoef([data[:,3],data[:,4]])
    print "MOD07press - CEM barometer R2:", b[1,0]*b[1,0]  
def corel_mod07temp_pm(time_range):
    f = open("correl/mod07temp_pm.csv", "w+")
    station_query = 'SELECT id, ST_AsText(location) as location FROM public.grdstation WHERE id >=1 and id <=5'
    mod_query="SELECT aqstime, ST_Value(rasref10, (st_GeomFromText('POINT({0} {1})', 4326)),true) as temp FROM apom.satresampmod07temperature where aqstime < '2014-09-10 00:00:00'::timestamp"
    pm_query="SELECT avg(pm1), avg(pm25), avg(pm10),avg(temp) FROM apom.grdpmhiscem_data WHERE (pm1 > 0 OR pm25 > 0 OR pm10 > 0) and aqstime > '{0}'::timestamp - interval '{1} minutes' and aqstime < '{2}'::timestamp + interval '{3} minutes' and stationid = {4} GROUP BY aqstime"
   
    con = None
    try:
        # doc danh sach tram
        con = psycopg2.connect(host='192.168.3.4', database='fimo', user='rasdaman', password='rasdaman') 
        cur = con.cursor(cursor_factory=psycopg2.extras.DictCursor)
        cur.execute(station_query)
        stations = cur.fetchall()
        for station in stations:
            stationid = station["id"]
            # bo qua tram Khanh Hoa
            # if stationid  == 3:
            #     continue
                
            location = loads(station["location"])
    
            cur.execute(mod_query.format(location.x, location.y))                                                                                                                                                                                                                  
            rows = cur.fetchall()
            
            for row in rows:
                temp_aqstime = datetime.datetime.strptime(str(row["aqstime"]), "%Y-%m-%d %H:%M:%S") + timedelta(hours=7)
                temp = row["temp"]
                if(temp!=None):
                    temp_value = temp
                else:
                    temp_value=0;

       
                if temp_value > 0:
                 
                    cur.execute(pm_query.format(temp_aqstime, time_range, temp_aqstime, time_range, stationid))
                    
                    pm_row = cur.fetchone()
                    if not(pm_row is None):
                        pm1 = pm_row[0]
                        pm25 = pm_row[1]
                        pm10 = pm_row[2]
                        cem_temp=pm_row[3]
                      
                        f.write("{0},{1},{2},{3},{4},{5},{6},{7},{8}\n".format(stationid, temp_aqstime.year, temp_aqstime.month, temp_aqstime, pm1, pm25, pm10,cem_temp,temp_value))
            
    except psycopg2.DatabaseError, e:
        print 'Error %s' % e    
        sys.exit(1)
    finally:
        if con:
            con.close()
        if f:
            f.close()
    data = np.loadtxt("correl/mod07temp_pm.csv", delimiter = ",", usecols=(4,5,6,7,8))
    print "Number of sample:",len(data)
    
    a= np.corrcoef([data[:,0],data[:,4]])
    print "MOD07temp - PM1 R2:", a[1,0]*a[1,0]
     
    b= np.corrcoef([data[:,1],data[:,4]])
    print "MOD07temp - PM25 R2:", b[1,0]*b[1,0]
     
    b= np.corrcoef([data[:,2],data[:,4]])
    print "MOD07temp - PM10 R2:", b[1,0]*b[1,0]
    
    b= np.corrcoef([data[:,3],data[:,4]])
    print "MOD07temp - CEM temp R2:", b[1,0]*b[1,0]

def corel_myd07press_pm(time_range):
    f = open("correl/myd07press_pm2.csv", "w+")
    station_query='SELECT id, ST_AsText(location) as location FROM public.grdstation WHERE id >=1 and id <=5'
    mod_query="SELECT aqstime, ST_Value(rasref10, (st_GeomFromText('POINT({0} {1})', 4326)),true) as pressure FROM apom.satresampmyd07pressure where aqstime < '2014-09-10 00:00:00'::timestamp"
    pm_query="SELECT avg(pm1), avg(pm25), avg(pm10),avg(barometer) FROM apom.grdpmhiscem_data WHERE (pm1 > 0 OR pm25 > 0 OR pm10 > 0) and aqstime > '{0}'::timestamp - interval '{1} minutes' and aqstime < '{2}'::timestamp + interval '{3} minutes' and stationid = {4} GROUP BY aqstime"
    con = None
    try:
        # doc danh sach tram
        con = psycopg2.connect(host='192.168.3.4', database='fimo', user='rasdaman', password='rasdaman') 
        cur = con.cursor(cursor_factory=psycopg2.extras.DictCursor)
        cur.execute(station_query)
        stations = cur.fetchall()
        for station in stations:
            stationid = station["id"]
            # bo qua tram Khanh Hoa
            # if stationid  == 3:
            #     continue
                
            location = loads(station["location"])
    
            cur.execute(mod_query.format(location.x, location.y))                                                                                                                                                                                                                  
            rows = cur.fetchall()
            
            for row in rows:
                press_aqstime = datetime.datetime.strptime(str(row["aqstime"]), "%Y-%m-%d %H:%M:%S") + timedelta(hours=7)
                pressure = row["pressure"]
                if(pressure!=None):
                    pressure_value = pressure*0.100000001490116
                else:
                    pressure_value=0;

       
                if pressure_value > 0:
                 
                    cur.execute(pm_query.format(press_aqstime, time_range, press_aqstime, time_range, stationid))
                    
                    pm_row = cur.fetchone()
                    if not(pm_row is None):
                        pm1 = pm_row[0]
                        pm25 = pm_row[1]
                        pm10 = pm_row[2]
                        barometer=pm_row[3]
                      
                        f.write("{0},{1},{2},{3},{4},{5},{6},{7},{8}\n".format(stationid, press_aqstime.year, press_aqstime.month, press_aqstime, pm1, pm25, pm10,barometer,pressure_value))
            
    except psycopg2.DatabaseError, e:
        print 'Error %s' % e    
        sys.exit(1)
    finally:
        if con:
            con.close()
        if f:
            f.close()
    data = np.loadtxt("correl/myd07press_pm.csv", delimiter = ",", usecols=(4,5,6,7,8))
    print "Number of sample:",len(data)
    
    a= np.corrcoef(data[:,0],data[:,4])
    print "MYD07press - PM1 R2:", a[1,0]*a[1,0]
     
    b= np.corrcoef(data[:,1],data[:,4])
    print "MYD07press - PM25 R2:", b[1,0]*b[1,0]
     
    b= np.corrcoef(data[:,2],data[:,4])
    print "MYD07press - PM10 R2:", b[1,0]*b[1,0]
    
    b= np.corrcoef(data[:,3],data[:,4])
    print "MYD07press - CEM barometer R2:", b[1,0]*b[1,0]  
            
def corel_myd07temp_pm(time_range):
    f = open("correl/myd07temp_pm.csv", "w+")
    station_query='SELECT id, ST_AsText(location) as location FROM public.grdstation WHERE id >=1 and id <=5'
    mod_query="SELECT aqstime, ST_Value(rasref10, (st_GeomFromText('POINT({0} {1})', 4326)),true) as temp FROM apom.satresampmyd07temperature where aqstime < '2014-09-10 00:00:00'::timestamp"
    pm_query="SELECT avg(pm1), avg(pm25), avg(pm10),avg(temp) FROM apom.grdpmhiscem_data WHERE (pm1 > 0 OR pm25 > 0 OR pm10 > 0) and aqstime > '{0}'::timestamp - interval '{1} minutes' and aqstime < '{2}'::timestamp + interval '{3} minutes' and stationid = {4} GROUP BY aqstime"
    con = None
    try:
        # doc danh sach tram
        con = psycopg2.connect(host='192.168.3.4', database='fimo', user='rasdaman', password='rasdaman') 
        cur = con.cursor(cursor_factory=psycopg2.extras.DictCursor)
        cur.execute(station_query)
        stations = cur.fetchall()
        for station in stations:
            stationid = station["id"]
            # bo qua tram Khanh Hoa
            # if stationid  == 3:
            #     continue
                
            location = loads(station["location"])
    
            cur.execute(mod_query.format(location.x, location.y))                                                                                                                                                                                                                  
            rows = cur.fetchall()
            
            for row in rows:
                temp_aqstime = datetime.datetime.strptime(str(row["aqstime"]), "%Y-%m-%d %H:%M:%S") + timedelta(hours=7)
                temp = row["temp"]
                if(temp!=None):
                    temp_value = temp
                else:
                    temp_value=0;

       
                if temp_value > 0:
                 
                    cur.execute(pm_query.format(temp_aqstime, time_range, temp_aqstime, time_range, stationid))
                    
                    pm_row = cur.fetchone()
                    if not(pm_row is None):
                        pm1 = pm_row[0]
                        pm25 = pm_row[1]
                        pm10 = pm_row[2]
                        cem_temp=pm_row[3]
                      
                        f.write("{0},{1},{2},{3},{4},{5},{6},{7},{8}\n".format(stationid, temp_aqstime.year, temp_aqstime.month, temp_aqstime, pm1, pm25, pm10,cem_temp,temp_value))
            
    except psycopg2.DatabaseError, e:
        print 'Error %s' % e    
        sys.exit(1)
    finally:
        if con:
            con.close()
        if f:
            f.close()
    data = np.loadtxt("correl/myd07temp_pm.csv", delimiter = ",", usecols=(4,5,6,7,8))
    print "Number of sample:",len(data)
    
    a= np.corrcoef([data[:,0],data[:,4]])
    print "MYD07temp - PM1 R2:", a[1,0]*a[1,0]
     
    b= np.corrcoef([data[:,1],data[:,4]])
    print "MYD07temp - PM25 R2:", b[1,0]*b[1,0]
     
    b= np.corrcoef([data[:,2],data[:,4]])
    print "MYD07temp - PM10 R2:", b[1,0]*b[1,0]
    
    b= np.corrcoef([data[:,3],data[:,4]])
    print "MYD07temp - CEM temp R2:", b[1,0]*b[1,0]
if __name__ == '__main__':
  
    #corel_mod04_aeronet(30)
    corel_myd04_aeronet(30)
    
    #corel_mod04_pm(30)
    #corel_myd04_pm(30)
    #corel_mod07press_pm(30)

    #corel_mod07temp_pm(30)
    #corel_myd07press_pm(30)
    #corel_myd07temp_pm(30)
